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I. INTRODUCTION 


Claims by many researchers indicate that a carefully 
designed controller could reduce fuel consumption by mini- 
mizing the propulsion losses which are caused by added drag 


due to steering of the ship. 


The goal of this thesis was aimed at developing апа 
demonstrating the utility of an improved steering control 
for the Mariner class ship. The immediate goal of this 
research was to develop design methodology for an adaptive 
autopilot that would provide effective steering control with 
associated cost savings for a full range of Seaway апа 


зна У conditions . 


To simulate the ship in the computer program, ship s 
nonlinear equations of motion were needed. Chapter 2 
addresses the Mariner class ship nonlinear equations of 


motion. 


Chapter 3 addresses the work on testing the ship simula- 
tion model and open loop ship's’ behaviors in calm water and 


in a seaway. 


The basic Nomoto models give an adequate description of 
ship steering dynamics for design. The Nomoto second and 
third order models were developed from the ship's’ linear 
equations of motion, Chapter 4 adresses the Mariner class 
Ship Nomoto model derivations by using mathematical methods 


and by using a function minimization subroutine. 


Chapter 5 shows an adequate cost function which repre- 
sents the added drag due to steering and includes deriva- 


tions for evaluating the weighting factor. Also Chapter 5 


LAE 


presents the assumptions and approaches needed to find the 


cost function that is used by many researchers. 


Ship dynamics change with operating conditions such as 
ship Speed, encounter angle and encounter frequency. 
Chapter 6 presents optimal controller parameters as a func- 


tion of different operating conditions. 


Chapter 7 addresses an approach to an adaptive controller 
utilizing information which is easy to measure on ship board 
such as ship speed, heading error and rudder angle. PHIS 
adaptive controller must be used to provide minimum added 


drag due to steering. 


Conclusions were drawn from simulation results I MIL 


presented in Chapter 8. This chapter also recommends some 
future studies, which can be done as extensions of this 
work. 
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II. NONLINEAR EQUATIONS OF MOTION 


Nonlinear equations of motion are suitable for predicting 
tight maneuvers and also suitable for computer programming. 


The nonlinear equations of motion used in this’ work have 


been developed by Abkowitz [Ref. 1, 2], and Strom Tejsen 
[Ref. 3], based on a Taylor series expansion of forces and 
moments. Terms higher than third order are not included in 


the equations because experience has shown that accuracy is 


not significantly improved by their inclusion. 


A result of symmetry about the xz-plane, X is an even 
Ewmetion of V, R, D, V and R so on, the crosscoupled terms 
in the equations involving odd powers of V, R and р аге 
zero, however, crosscouple terms which involve even powers 
of V,R and D are nonzero. In contrast to X, the expressions 
for Y and N are odd functions of V,R,D,V and R: пакты 
only the coefficients of the terms in the expansion with odd 
powers are nonzero; those with even powers are zero. For 
Some reasons, Х 15 neither an odd nor an even function of U 


but rather its expansion includes all powers of DU. 


Equations X,Y and N are functions of OR IJ V R and D. 
Taylor series expansions of X,Y and N including terms up to 


the third order are as follows: 


(к-сі - £(U,V,R,D) (еап 2.1) 
(m - Y,)*V * (m*X,- Y; )"R = f(U,V,R,D) (eqn 2.2) 
(m*X, - Nea + (LS N R = £(U,V,R,D) (eqn 2.3) 
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Where: 


ус алуа e Хау: d^X/dU*dV 
= = 3 2; 
Ке @?X/dU? Xu ААО ае 
Е(П,У,К,р) = 0.5% #0? + (1/6)#Х а УИ (eqn 2.4) 
. 0,555, УІ: (ШЕ a 
+ 0.5%X Алуа + 0.58Х ла ВРЗОО 
- (X oa m)*V*R + ов + Хум КОШ 
+ Хос V DDU + Xoga RKD DU + Nx DU + me 
+ 0.5%*Xo,*D? + 0.5%X gq D? *DU * X4 *R*D 
и Serene © о 
f(U,V,R,D) - а "В Үз Бат ты» (egy 255) 


+ Y¥° + O.5*Y na VR? + ООо Ер 
* Y4uAV*DU + 0.58 ууъ УРОШ? + (Ма- пеша) «Е 
* (1/6)*Y, n *D? * Yea ЫП 2-5 О. 


e 0.5*Y aR DR t You DADU + © Би 
| о УКР + ( ey, 6 даде “ү 3 4 0. 5*Ya pp" R*D 2 


£(U,V,R,D) s Nt, *DU « NA, *DU* « Ny*V « Nj*D — (eqn 2.6) 
+ ОМ. DADU S (No -m*X,U1)*R 
+ МУ" + ОБАМА ауы М иво УКР + № 
к 0.55 РОМ? + 0.52Мрра 2982 + Ny y*D*DU 
+ NgužR*DU + 0.5% ВРО? + (1/6)2Мррр 
*(1/6)*Na R + б.р Ори вв. 
+ (116) Noy V ИВО ОТУ ТУ. ODENSES V Di 


All of the derivative coefficients of the equations are 
evaluated on the basis of experimental data obtained from 


captive model tests, and given in Table I, II and III. 


The Y force and N moment induced by the rotation of a 


single propeller or by unirotating multiple propellers at 
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V=D=0.0 are identified as Y and № these are likely to be 
Speed dependent. To see the rudder effects in calm water or 


Sea state propeller effects were ignored. 


Finally, from the X , Y and N equations the ship's surge, 


sway and yaw equations can be written as follows. 


£(U, V, R,D) 
ЕС o Е leagan 2-7) 


(m - Xi) 


u 


(T,- N, )*£(U,V,R,D) - (m*Xc- Yi )*£(U, V,R,D) 


<. 
ii 


(eqn 2.8) 


м 


2) G - Y.) 


(м - u s; (m*X. - Na )*£(U, V, R, D) 
D s O—— — (еап 2.9) 


ви г). №) - (m*X. - Ya) 
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These surge, sway and yaw equations can be rewritten in the 


form: 
dU/dt = g(t, U(t), VE) R(T PSD r 


dV/dt Blt, ОЕ), ЕЛ, ВЕР (eqn 2.10) 


dR/dt e(t; U(t); (Е), RE) ИШ 


Where U(t), V(t), R(t) and D(t) are the instantaneous values 
of U, V; R and D at any eime En 


Equation 2.10 is a set of three first-order differential 
equations for which approximate numerical solutions are 
readily obtained on a digital computer, The key to the 
numerical solution is that values of U, V and R at time 
t+DELT are obtained from knowledge of the values of U, V, К 
and D at time t using a simple first-order expansion; that 
15, 

U(t+DELT) = U(t) + DELT*U(t) 


V(t+DELT) = V(t) + DELT*V(t) (eqn 2.11) 


R(t+DELT) = R(t) + DELT*R(t) 


This method is found to give adequate accuracy for the 
present type of differential equations because of the fact 
that the accelerations U, V and R vary but slowly with time, 
due to the large mass and inertia of a ship compared to the 
relatively small forces and moments produced by its control 
surface. Any desired accuracy of the solutions сап Бе 
obtained with a computer by simply using smaller time inter- 
vals DELT. This procedure was used for all computer programs 


which were developed for this thesis. 
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TABLE 1 


Aecesccmemi of the Coefficients’ in 


Variable 
U 
DU 
DU? 
DU? 
V? 
R? 
D? 
V**DU 
R**DU 
D**DU 
V*R 
Мер 
R*D 
V*R*DU 
V*D*DU 
R*D*DU 


Coefficient 


Value 
По 

e 

ВВ Хх 

ШЕ. 

su Gy 

о Xo 

Боа 


к 


О Оо о О 


OX ans 
0.5*X 9.0 
Худ“ т 

Х 49 

XRD 
Хурал 
Хуох 


X aou 
n 


the X Equation 


ОЕ бое ЕЕ стеле 
0 177 
20.0253 
0.00948 
-0.00217 
-0.189 
(0001037 9 
- 0102 


O 168 
0. 0196 


Ка ел са тег аге попа!тепгтопа|1г7еада on the basis of 


RHO,L,T and S 


No entry in these columns means the coefficient was ignored. 


I; 


Variable 
V 
R 


*R? 
“р 2 
*DU 
*DU? 


мј с< < < < < < 


ЕЗ 


R*D? 
R*DU 
R*DU? 


D? 
ру? 
D*R? 
D*DU 
D*DU? 
УВ 


DU 
ЕШ 


ТАВБЕ ФИ 


Assessment of the Coefficients in the Y Equation 


Coefficient 
nom 


coo и 
0. К 
0.5- Y 
Yuu 
Оса пе 
(Yer m ) 
(1/6 )*Yo oe 
Па oe 
0.5 99 


VOD 


Yom 
0.5*Y, iu 
Yo 
(1/6)#Yooo 
0 е 95 
ОУ 
You 
0.5*Y uu 
Хчко 

ya 

Ша 


Y A 


Value of Coefficient 
E327 
-0 001 
-0.244 
-1.702 
0 
-0.0008 


50: 105 


9228 


0.0586 
= 0.00975 
075 

0 
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ШАРЕ II 


Assessment of the Coefficients in the N Equation 


Variable Coefficient Value of Coefficient 
V (m*X.- N.) -0.00478 
R (1,- %) 0.0175 
V Ny -0.0555 
y? Е 0.345 
V*R? ONSE. 0 
ухо: 0.5*N o 0.00264 
ухру Noy _ 
VDU? ОБЕМ лула, з 
R (Ny- m*Xg*U1) -0.0349 
R? (1/6) “Мах 0 
R*y? OS Nou БЕ 
R*D? м 0 
R*DU Neu ú 
R*DU? 0.5*Na uu _ 

D No 50.0298 
D (1/6)*N po 0.00482 
p*y? Он, 0.1032 
D*R? 0.5“ ар 0 
D*DU Nou 0 
D*DU? 0.5 *N А а 
V*R*D Nee 0 

_ № 0.00059 
DU Мо, 0 
DU? № B 
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А. CALM WATER CASE 


Using a ship's nonlinear equations of motion, and Mariner 
class ship coefficients, a simulation program THESIS FORTRAN 
was developed and run to observe U, V and R.The computer 
program THESIS FORTRAN is shown in Appendix A. 


First run D-0, YAWC=0, Ul-15 knots were applied and it 
was seen that the ship stays on its initial course апа 
speed. U=15 knots, V=R=0. 


The program was rerun a few times, changing the rudder 
angle to 2.8 degrees to both sides (port and starboard) and 
it was observed that increasing the rudder angle changes the 
ship's course and also U decreased, the absolute values of V 


and R increased. After a few hundred seconds U, V and R 


reached steady values independently. These steady values 
depend on rudder angle, large rudder angles decrease U and 
increase V and R. For a rudder angle of one degree and 
constant speed (Ul) of 15 knots, the first 200 seconds of 
time response of U, V and R are shown in figure 3.1 as an 
example. 


В. БЕА STATE CASE 


To observe the behavior of the ship in a sea state, 
disturbance forces and moments are needed that depend on sea 
state, ship speed, encounter angle and encounter frequency. 
Also in sea steate hydrodynamic parameters are changed, i.e, 
the added mass and added inertia are functions of encounter 


frequency and sea state. 
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The disturbance forces, moments, and added mass, added 
inertia terms were found by running the sea state program 
that is presented in [Ref. 4]. Data for a Mariner class 
ship and chosen sea conditions that were used in the sea 


state program are shown in APPENDIX B. 


For observation purposes 'FX' (disturbance force in 
Surge),  'FY' (disturbance force in sway) and 'MZ' (distur- 
bance moment in yaw) were added into the surge, sway and yaw 


equations that were used in the simulation program. 


For regular seas the program has been run a few 
times.Every time different FX, FY, MZ and coefficients were 
used to represent different ship characteristics such as 
ship speed, loading * etc. and enviromentel conditions such 
as sea state, encounter angle, encounter frequency. Results 
Show that U, V and R are sine waves with amplitude and phase 


depending on ship and environmental conditions. 


Time response of U, V and R in 200 seconds are presented 
in Figure ( 3.2), for ship speed 15 knots, encounter angle 
030.0 degree, encounter frequency 0.60 radian/second and sea 
State 6. 


*During this research, displacements (Molded) up to 32 
ft were chosen as a loading condition for the Mariner. 
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Figure 3.2 Time Response of U,V and R in Regular Seas. 
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IV. NOMOTO MODEL OF TAERA RINER 
To find a model which can be used in computer simulation 
the Mariner's linear equations of motion and its hydrody- 
namic coefficients were used and third and second order 
Nomoto transfer functions derived. Values used were for ship 
speed of 15 knots. 


Mariner s linear equations of motion are 


(m - X. )%U = X. *DU (eqn 4.1) 


“үу = (У. - 


(тп т Y, )*V - Y R 


y m*X¢)*R $ Coe m*U1)*R 


el UN: 


5 )*R 2. (М, - mke UD R = (Ny - m*Xe )*V ve У 


А. MATHEMATICAL APPROACH 
Proceeding to the third order Nomoto equation: 


YAW(s ) КОСИБ, 


D(s) Ss (Pls +I) P Пати) 


Deriving the third order Nomoto transfer function from 


the sway and yaw equations, we show them in matrix notation 


as follows. 


V = 03572 о V БГИ 


R -.0003 -.10 -. 
Y 10 | |R 003 | 
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X = A*X + B*U then, X(s) - (s*I-A)!'!*B*U (assuming initial 


conditions are zero.) 


V .2101^5 + .0534 





ЕЕ) (116.93%s+1) 


R -.0038%5 - .0002 


R(s) = s*YAW(s) 


YAW(s) `` -0.189*(18.34*$+1) 
So, then — n = 


D(s) s*(7.67*s*1)*(116.93*s«1) 


чине K=O. 189, 2=18.34, Pl=7.67 and ?2-116.93 


Proceeding to the second order Nomoto equation: 


YAW(s ) K 


D(s) S*(Pl*s + 1) 


Deriving the second order Nomoto transfer function from 


the yaw equation only, the result is K-0.03 and P1-10. 


в. COMPUTER APPROACH 


Ме used a function minimization Subroutine to obtain 


parameters of the transfer functions. Figure 4.1 shows the 
scheme used to obtain the third and second order Nomoto 
transfer functions. The results are given in Table IV. 
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Figure 4.1 Determination of Third and Second Order 
Nomoto Models. 
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TABLE IV 


The Nomoto Model Parameters for Mariner 


Third Order Second Order 
Е = 0.189 K =0.0298 

zm 18367 D-90599 
PEDE 7.6739 
БРА 1100929 


As Seen the answers obtained by function minimization 


agree closely with the analytic solutions. 
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V. DERIVATION OF A COST FUNCTION FORTHE SDD I AE P P Tp 

In recent years many researchers have studied the problem 
of optimizing an automatic ship steering controller for 
minimum fuel consumption. It 1s well known that additional 
drag is introduced by steering and that both the rudder 


motion and the yawing motion contribute to this added drag. 


When deriving a cost function we are required to find one 
cost function that must be convenient for ship board use.The 


cost function that is commonly used in recent years is 


t 


o 


J = | CEAMDA“YAWE ~ ка Бр Јосе (eqn о =) 


Where YAWE = Yaw error 
LAMBDA = Weighting factor 


Derivation of this cost function from the surge equation has 


been well explained in [Ref. 5] for the SL/7 class ship by 
RES Ел 


To derive the cost function for ithe Marinei REID's 
approach is taken as a reference and his assumptions are 
used. To show how the cost function for the Mariner was 
derived, REID's work is presented here step by step with 


derivations of the Mariner's cost function. 


The Surge Equations: 


SE T 
(m - Xa) U ОЗИ 0.5 Di (eqn 5.2) 
+ (Xyg + m)*V*R + Хр 
Where: К - -0. ООО ОБЕ 


Хр Propeller thru rt 
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МАКТМЕК: 


(m - X4,)*U » 0.5*X,,*V? « (Xyg * m)*V*R (eqn 5.3) 
* 0.5*Xg "D? * (0.5*Xq. * m*XgQ)*R? 
О 5х0? 


(ШКО) X DUE 


+ 


+ 


МА 


It is seen that there are some terms in the Mariner's 
Surge equation which they are not included in the SL/7's 
surge equation. Assuming steady state situations since, U-O0. 


The instantaneous surge relevant to steering is 


ЕК - 0 55Х U 0.5*X. 05 + (Xj. * m)*V*R (eqn 5.4) 


SEI: 


MARINER: 


Е ОЕ о + (Хо + m)*V*R (eqn 5.5) 
+ (0,584 + m#X )*R? + Xyp*V"D + X.#DU 
+ 0.5*X4,*DU? *(1/6)*X ОО? 


From the instantaneous surge equation relevant to 


Steering of the SL/7 Reid came up with the following cost 
СО: 


€ 


о 
J = (1/2T)* / (LAMBDA"*V*R + ETA*V? + D?)*dt (eqn 5.6) 


о 
Where LAMBDA"- (m * X 48) /0. 5*X gg 


р ОО Бах 


29 


and һе found that the value of the (ETA*V?) term is very 
small so he neglected the (ETA^V^) term, then the cost func- 
tion for thes5p 15 


J = 0.5*/(LAMBDA"*V*R + D?)*dt (eqn 5.7) 





Using the same approach for the Mariner, the following 


cost function was derived. 





J = (1/2T)*/(AL*DU * A2*DU? + A3*DU? + A4*V2 — (eqn 5.8) 


- AS*R? + р? « AZ*V*R - A8*V*D)*dt 


Where № = К ORE А2 = 02.59 S ОР ти 
АЗ = (1/6)*X\../0.5%X5, АА = 0. | 
дБ = О ао wen 
А7 = (Ха % п)/0.55Хзу А8 » Xyg/0.5* Xs 


A5 and A8 are always negative numbers, since then A5 and A8 
have a minus sign in the equation. For the calm water case 
and when Ul=15 knots, D-2.6 degrees after 2000 seconds, 
values of every term in the Mariner's cost function are 


given below to give an idea about assumptions. 


Al*DU = -0.001939 A2*DU? = -0.00000111 
A3*DU?- -0.00000000039 А4Х 2 = 0.000003253 
А5*В? = -0.0002884 АДЕМ =s 0050001923 
р? = 0.002059 A8*V¥D = -0.00002609 


As seen from the above (A2*DU? ), (A32DU ), ЦА) аа 


(A8*"V*D) terms are very small compared to others, so they 
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may Бе neglected. Also to measure the DU оп shipboard is 
very hard although it may someday be measured by new satel- 
ПЕЦЕ Кае тес со уе must mot include it in the cost func- 
tion. After these assumptions the cost function for the 


Mariner is 


B 


o 
ew >< АВ + р —* AZ*V*R)*dt (еап 5.9) 
0 

The only difference between Equation (5.9) and Equation 
(ти тле (25 |) term that 15 пос included in the SL/7's 
cost function. To see the effect of the (A5*R?) term on the 
Eus uU—Atlon' the Mariner s cost function was evaluated with 
and without the (A5*R?) term for the calm water case and 


Ul-15 knots, D-2.6 degree after 2000 seconds, results are: 


J 
with (AS*R?) 0.0025399826296 
without (A5*R?) 0.0022515066462 


— = аә $ -. -. $ se $ $5 $ - кы. -. |= =— -. BF -. «те «е «е «е BSF «е = = = = = = = = = = = = = = 


There is no big difference between these two J values, and 
to make the derivations similar to Reid's derivations the 
Егет чоп t be included in the Mariner's cost func- 
tion but as it is known that for the Mariner the (A5*R*) 
term is as big as the (A7*V*R) term, it would be better to 
eons der It im the cost function. After all of the above 
Steps the cost function of the Mariner may be written as in 
Eq@ation (5.7). 


JL 


V and R are hard to measure on shipboard, Burts.) NUM 


be defined as 
V*R = OP*YAW? 
Where R = YAW = YAWE*w 


OP = Distance from the ship pivot point to 
the origin. s.0.3*L 


w - Natural frequency of the ship's steering 


system closed loop. (w - 0.05 rad/sec was initially used.) 
Finally the cost function for the Hariner is 


+ 


5 
Ј = 0.5% |(ЏАМВРАХУАМЕ2 * D?)*dt (еап 5.10) 


° 
Where LAMBDA = (m + Xua )%*OP%*w2/0.5%Xyy 


Since, X depends on ship speed, for different ship speeds 
values of LAMBDA were calculated and are presented in Table 
(У). These LAMBDA values were calculated by assuming the 
natural frequency of the ship's steering system closed loop 
is equal to 0.05 radian/ second. How important the accuracy 
of the LAMBDA value is with respect to finding the optimal 


control parameters will be observed in Chapter 6. 
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TABLE И 


Values of ELA Factor for Different Speeds 
The Hariner Class Ship 


же = = = шш «< = = == = = = = = = = = = <= => = = = = = «кы. <= <= = = = = = == ж. <= = <= <= <= = <= = = “ы “ы “ы «ы == 


ship speed | LAMBDA 
(Knots) | 
p" i j 65 
E a a 
N СЕ 


33 


VI. CONTROLLER DESIGN FOR THE MARINER, USING FORTRAN PROGRAM 
We coupled a function minimization subroutine to the 
cascade compensater that is coupled with Mariner's equations 
of motion and used the subroutine to adjust the controller 
parameters to minimize the cost function which is derived in 
Chapter V, and evaluate the minimum cost.The Fortran program 
to calculate the optimal parameters ыы Appendix C. 
Compensater ‘A' and 'B' are used as controllers, their 


structures are shown in Figure ( 6.1 ). 


Кх(5х21 + 1) Kx(SxZl + 1) 


(SxP1 + 1) (SxP1 t 1)x(SxP2 t1) 


COMPENSATOR- A COMPENSATOR ` 8 





Figure 6.1 Various Structure Controllers: 


Aw, CALMAWATER CASE 


For calm water, for a given 1 degree yaw command, using 
the computer method we optimized controllers 'A' and 'B' and 


the results are shown in Tables Vi V 
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ТАВЬЕ УТ 


Optimal Controller Parameters I 


Similatiom Results - Steady State 600 seconds 


For ship speed 10 Knots the optimal 
parameters of various controllers 
and the cost 


CONTR K za Pl P2 COST 

А 0.77 60.07 20.66 = 0.0494896 

B 0.74 60.07 2912 0.902 0.0514841 
TABLE VII 


Optimal Controller Parameters ТТ 


Simulation Results - Steady State 600 seconds 


For ship speed 15 Knots the optimal 
parameters of various controllers 
ana the cost 


CONTR K p Pl P2 COST 
А 035 51.46 18.08 ----- 0.0186974 
В 0.50 51.58 17.81 0.890 0.01957901 
TABLE VIII 


Optimal Controller Parameters III 


Simulation Results - Steady State 600 seconds 


For ship speed 20 Knots the optimal 
parameters of various controllers 
and the cost 


CONTR K A Pl P2 SOSH 
A 0.40 44.87 REO E. 0.0094217 
B 0.39 Дт 15.84 0.880 0.00991801 
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These results will be references forthe contro ЕЕ 
design for sea state operation. We observe that increasing 
the speed gives us smaller controller parameters, also this 
behavior can be seen from the о „| ке пе ЕЕ е 


presented in [Ref. 6]. 


B. REGULAR SEAS CASE 


The ship in regular seas is affected. by sea wave distur- 
bance forces and moments. These are functions of sea state 
and encounter frequency. Also the added mass and the added 
inertia terms are functions of sea state апа encounter 
frequency, and the encounter frequency depends оп the 
encounter angle and ship speed. All of these variables must 
be considered when calculating the optimal parameters of the 


controller: 


Using the computer method controllers 'A' and 'B' were 
optimized for a few different cases and the results are 
shown in Tables IX, X, XI, XII, XIII; XIV, KV: T апа ый 


TABLE IX 
Optimal Controller Parameters IV 
Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 


For ship speed 15 Knots, encounter MS 030.0 па ше, 
encounter frequency 0.50 rad./sec. and sea state 6. 


CONTR K zi Pl P2 COST 
A 0.358 66.6 24.61 ----- 0.001252939 
B 0.35 44.68 06.58 9.720 0.0010086581 
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TABLE 8 
Optimal Controller Parameters V 
Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 


For ship speed 15 Knots, encounter ane c 060.0 SEES 
encounter frequency 2.50 rad./sec. and sea state 9. 


CONTR K zii P1 P2 COST 
A 0.60 41.58 ШЕЕ —-— = 0.000020747 
В 0.70 30.49 ОЗИ 2.940 0.000016038 
TABLE XI 


Optimal controller Parameters VI 
Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 


For ship speed 15 Knots, encounter ana 090.0 asss ss 
encounter Irequency 0.50 rad./sec. and sea state 7. 


CONTR K Zl Pl P2 COST 
А 0.54 50.65 09.35 ----- 0 0542239898 
B e ed sS se sS s< IT DID NOT CONVERGE Е 
TABLE XII 


Optimal Controller Parameters VII 


Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 


For ship =: 15 Knots, encounter апе 2070 ее 


encounter frequency 0.75 rad./sec. and sea state 
CONTR K zu Pl P2 COST 
А 0,52 41.09 08.95 ----- 0.018345 
B s sss IT DID NOT CONVERGE > 
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TABLE XIII 
Optimal Controller Parameters VIII 


Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 


For ship ` 15 Knots, encounter an 15070 ERTES 


encounter frequency 1.50 rad./sec. and sea state 
CONTR К Zi P1 P2 COST 
А 0.645 42.40 07:70 == 0.0008188 
В еек TT DID NOT CONVERGE ккеккиеккиикии 
TABLE XIV 


Optimal Controller Parameters IX 


Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 


For ship еве 15 Knots, encounter ` Е 090.0 s tee 


encounter reguency 0.50 гаа. / ѕес. ап sea starce 
CONTR K 21 Pl P2 COST 
А 033 27.45 03704 2... 00/2869 79/59 
P VoU IT DID NOT CONVERGE уеетеуелеугзе еле сус 
ТАВЬЕ Ху 


Optimal Controller Parameters X 
Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 


For ship speed 15 Knots, encounter ame he 090.0 degree, 
encounter frequency 0.50 rad./sec. and sea state 6. 


CONTR K zu Pl P2 COST 
А 0256 SION 11.87  ----- 0.019029424 
В о IT DID NOT CONVERGE WoW а а 
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TABLE XVI 


Optimal Controller Parameters XI 


Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 


Bor ship S 10 Knots, encounter ШЕШЕ 060.0 о 


encounter frequency 2.50 rad./sec. and sea state 
CONTR K Zl P1 P2 COST 
А 0.80 45.00 10.66 |  ----- 0.000045194 
B 0.89 36.27 03522 03 22 0.000036077 
TABLE XVII 


Optimal Controller Parameters XII 


Simulation Results - Steady State 600 seconds 
optimal parameters of various controllers and the cost. 


For ship Eco 20 Knots, encounter е 060.0 degree, 


encounter frequency 2.50 rad./sec. and sea state 
CONTR K Zl PI Ва COST 
А 0.47 34.67 10.35 |  ----- 0.000014376 
B OST 23 Ще 03:723 03205 0.000011085 
ÀS Seen from the above results, for all cases the 


compensators have characteristics of a lead network. 


The effects of sea state on the controller parameters 
can be seen by comparing the Tables ХУ, XI and XIV AS Seen 
from the tables an increase in sea state causes an increase 
in the cost value and a decrease in the controller parame- 
ters as expected, because heavy sea state brings high 
disturbance forces and moments and they cause heavy yawing 


motions. To answer this, the controller time constants 
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decrease and new parameters introduce more rudder motion, so 


an increase in rudder motions and heading error causes 


increase in the cost. For ship speed 15 knots, encounter 
angle 90.0 degree, encounter frequency 0.5 radian/second, 
sea state 6 and sea state 8, rudder angles and heading 


errors were plotted in 200 seconds and they are shown in 
Figure ( 6.2 T mand (TEI 
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Figure 6.2 Rudder Angles in Degrees for Sea State 6 and 8. 


40 


'E) 


GRE 


Y 
4 


YAWE (DI 





PNE 001205287 32983609 410 451 492 533 574 
TIME (SEC.) 


Figure 6.3 Heading Errors in Degrees for Sea State 6 and 8. 


Пет зе seen from Figure ( 6.2 ) and ( 6.3 ), for sea 
State 8 rudder and heading error values are bigger than the 
sea state 6 rudder and heading error values, but it was 
noticed that the periods are almost the same for both sea 


Seaces. 
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To have a better understanding about the effects of sea 
state on optimal controller parameters the cost value versus 
parameter values were plotted for ship speed 15 knots, 
encounter angle 90.0 degree, encounter frequency 0.50 
radian/second and sea state 6 between 8 and they are shown 
in Figure ( 6.4 ), ( 65), (66) апт ШШШ n n 





Figure 6.4 The Cost vs. K. 
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11 


Dum 


12 


13 


Curve "А and В аге projection of main curve 
ӨП ООЛ апа JPI surfaces. 





Figure 6.7 ieee eCasteyvs.w2) and Pl. 


As seen from the figures the optimal parameter values 
have a linear behavior when changing the sea state. This 
property would be useful to create look up tables for param- 


eter values. The use of look up tables will be discussed in 
Chapter 7. 


To see the motion of the ship, the ship's equations are 
solved in ship coordinates (x , y) and then transformed to 


Space coordinates (x , y), Figure НИКЛА поне спейоттепса 


tion of space axes and moving axes. 
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POSITION OF CENTER OF GRAVITY 
OF SHIP AT TIME to 


0 SPACE AXIS 











TRANSVERSE 
POSITION 
COORDINATE 2 


YOG Lf 
и А 
\ Ú (АКО 
N и \ 
PN 2220. | 
ОС cU N 
LONGITUDINAL ~~ 7 \ 
COORDIN TE N = 
COORDINA m 7 
au VELOCITY VECTOR 
TANGENT TO 
INSTANTANEOUS 
SHIPS PATH 
Yo 
+y 
Figure 6.8 Orientation of Space Axes and Moving Axes. 


The transformation from ship to space coordinates 1$ 


defined as 


Xog - U*cos(YAW) - V*sin(YAW) (eqn 6.1) 
U*sin(YAW) - V*cos(YAW) 


Yog 
Where Xog, Yog - Coordinates of the center of mass of 


the ship relative to coordinate system 


fixed with respect to the surface of the 


earth 
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ВЕ о was added imto the fortran program and 
the ship's motion were observed for 600 seconds in calm 
water and in regular seas for 30.0 degree encounter angle, 
sea state 6, encounter frequency 0.60 radian/second and ship 
speed 15 knots. Both compensator "А! апа "В" results are 
shonn in Figure ( 6.9 ), and after 600 seconds ship's coor- 


dinates are given in Table XVIII. 


TABLE XVIII 
Госагтон ОС the Ship 


Simulation Results - Steady State 600 seconds 


CASE Xog (ft) Yog (ft) 
Calm water 
No rudder 15200.994873 0.00 


Regular seas 
Comp.'A' 15 T99- 0237114 57:50 


Regular seas 
Comp. B- №5199. 06-219 2097116 


For calm water and regular seas compensator 'B' gave 
better results than compensator 'A' did, but for some cases 
compensator 'B' did not converge and the difference between 
results is not significant. The comparisons were made as to 
which type of compensator brings the ship nearest to final 


location at the end of the 600 seconds. 
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ҮОС (ЕТ) 





0 2000 4000 6000 8000 10000 12000 


14000 16000 


XOG (FT) 
Curve a ------ Using compensator ‘A’ 
Curves ТЕЕ Using compensator B 
Curve c --..-- The ship is in calm water 
(D=0 and YAWC=0) 
Бірптге 679 Simulation Results-Steady State 600 sec. 
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То see the effects of the weighting factor ( LAMBDA ), 
on the compensator parameters, different LAMBDA values, were 


used with compensator ‘'A' and the compensator parameters 


were computed for ship speed 15 knots, calm water and 
regular seas with encounter angle 90.0 degree, encounter 
frequency 0.5 radian/second and sea state 8. Results are 


shown in Table ( XIX) and Table ( XX ). 


TABLE XIX 
Test for LAMBDA Value I 


Simulation Results - Steady State 600 seconds 


for different weighting factor values, 
calm water, ship speed 15 knots, 
optimal parameters of AÀ type controller. 


LAMBDA K Zu РІ 
10.00 15077 29.59 о 
9.00 1.00 5 11.90 
8.00 0.94 32.91 1225 
7.00 0.88 34.84 13.14 
6.00 0.80 37.61 14.02 
Saal 0.76 39.38 14.62 
4.91 0.72 41.26 15.19 
me 0.68 S ПБ 8 
3.91 0.63 45.26 16.36 
Веда 0.58 48.06 17.10 
2.91 0.53 ОШО; 17.93 
2.41 0.48 2521 18.99 
1.91 г) 2% 59.79 PONES 
В 0.35 64.79 Т 
0.91 Doo 7 76.89 2 
RET 0.18 89.61 28 20 


&9 


TABLE XX 
Test for LAMBDA Value II 


Simulation Results - Steady State 600 seconds 


for different weighting factor values, 

regular seas, ship speed 15 knots 

encounter angle 90.0 degree, sea state 8, 

and encounter frequen 0.50 ra Есен 
У 


optimal parameters of pe controller. 

LAMBDA K Le E 
5.4] 0.71 24.91 08.41 
4.91 0.68 2556 08.37 
4.41 0.64 25589 08 -32 
3:91 0.60 26.34 (06522 
3.4] 0.56 26.85 08.11 
_ 21 0.55 27.45 08.04 
2.41 0.48 29528 08.04 
bcp 0.44 29 6 08.16 
1.41 0.40 3] 21 08.35 
0.91 0.34 33 а 08.71 
0.41 0.26 39.00 09 27 


Using above compensator values did not make а signifi- 
cant change on the ship location after 600 seconds simula- 
tion, so it can be said that the accuracy of the weighting 


factor 1S not very important. 


A few Simulation runs were performed by changing the 
optimal compensator parameters a small amount and the cost 
curve was plotted. Figure ( 6.10) shows the cost curves for 


three different K values versus Zl and Pl when the ship is 


in calm water, ship speed 15 knots and 1l degree course 
change, and Figure ( 6.11 ) shows the cost curves for ship 
Speed 15 knots, encounter angle 60.0 degree, encounter 


frequency 2.5 тайзап/зесопа, апа ве © 
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Changing Around the Optimal Values 
for the Sea State Case. 
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Ihe results showed that the cost curve is flat around 
the minimum. The same conclusion is made in [Ref. 7 and 8] 
fo the SL//. Using controller parameters on the flat 
Pemeron ot the cost curve, but not at the optimum did not 
meee a significant differences in simulation results, after 
600 seconds of cruising there was little change in the final 
location of the ship. This property of the cost surface may 
be useful in reducing the time required to find a pratical 


ШҮП for the cost function. 


The ship speed effects on compensator parameters for 
regular seas can be seen by comparing Tables X, XVI and XVII. 
Observe that for both cases increasing the ship speed made 
the parameter values and the cost value decrease. The same 
observation can be made by comparing the results that are 
found for the SL/7 container ship for different speeds іп 
[Ref. 9]. These results strongly indicate that the dynamics 
of the ship determines the optimum structure for the 


Controller. 


To see the stability of the system, the ship's third 
order Nomoto model was cascaded with the compensator and the 
open loop system BODE plot was drawn, in every case the 
system is stable. The phase margin varies between 40 
degrees and 70 degrees and the zero cross over of the magni- 
се ецшуе 15 аљоштпа 0.04 radian/second. ТЕ was observed 
that small changes in the compensator parameters do not 
affect stability. With in large limits of parameter varia- 
tion the system is always stable. Бот regular seas, ship 
speed 15 knots, encounter angle 30.0 degrees, encounter 
frequency 0.6 radian/second and sea state 6, using compen- 
sator 'A' and 'B', the structure of the system is presented 
in Figure ( 6.12 ) and the system open loop BODE plot and 
NICHOLS plot are shown in Figures ( 6.13 ) and ( 6.14 ). 
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F'ipguret6.12 Open Loop Steering Model. 
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Changing the environmentel conditions changes the 
optimal controller parameters. A few simulations were run 
changing the environmental conditions only but the optimal 
controller parameters were kept unchanged The behavior of 
rudder and heading error were observed. Three cases were 


defined depending on the enviromental conditions. 


Case I :Encounter angle 90 degree, encounter frequency 
0.5 radian/second, Sea state 6 and ship speed 15 knots. 


Optimal controller parameters for Case I are 
К = 0,56 Z9. 39 7]. Pl = 11.67 


Case II :Encounter angle 90 degree, encounter frequency 
0.5 radian/second, Sea State 8 and ship speed 15 knots. 


Optimal controller parameters for Case II are 
К = 0.53 21 = 27.45 Pl = 08.04 


Case III:Encounter angle 30 degree,encounter frequency 
0.5 radian/second, sea state 6 and ship speed 15 knots. 


Optimal controller parameters for Case III are 
K = 0.358 Zl = 66.60 DIM M NO 


Figure ( 6.15 ) and figure ( 6.16 ) show rudder and 
heading error when the ship is in Case I condition using 


Case I and Case II parameters. 


Figure ( 6.17 ) and figure ( 6.18 ) show rudder and 
heading error when the ship is in Case II condition using 


Case II and Case I parameters. 


Figure ( 6.19 ) and figure ( 6.20 ) show rudder and 
heading error when the ship is in Case III condition using 


Case III and Case I parameters. 


нише ( 6.21 штапа tipure ( 6.22 ) show rudder and 
heading error when the ship is in Case I condition using 


Case I and Case III parameters. 
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Figure 6.16 Heading Error in Case I with Case I 
and Case II Parameters. 
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Figure 6.17 Rudder Motion in Case II with Case II 
and Case I Parameters. 
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Figure 6.18 Heading Error in Case II with Case II 
and Case I Parameters. 
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Figure 6.19 Rudder Motion in базе TII with Case III 
and Case I Parameters. 
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Figure 6.20 RESCUE n се III with Case III 
as arameters. 
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Figure 6.21 Rudder Motion in Case I with Case I 
and Case III Parameters. 


61 


LEGEND 
USING CASE T PARAI GC A BI 
о USING CASE nee he) 


0.3... 0.5 0. 


0.1 


ы 
= 
oe 
gen. 
5 
z 
== 
= 
2 


О 


| і | i 
82 123 164 205 246 287 328 369 410 451 492 533 574 
TIME (SEC.) 





Figure 6.22 Heading Error in Case I with Case I 
and Case III Parameters. 


As seen from Figures ( 6.15 ), (6.16 ), (6.17 ) and 
( 6.18 ) using Case I optimal parameters in Case II or Case 
II optimal parameters in Case I did not make а big differ- 
ence in rudder motion and heading error, except in the tran- 


sient response part. 


Figures ( 6.19 ) and ( 6.20 ) show that to use Сс е 
optimal parameters when ship is in Case III is not proper 


because those parameters increase rudder and heading error. 


As seen from Figure ( 6.21 ), using Case III parameters 
in Case I provided better rudder motion after the transient 
response part. Cost values were calculated and it was 
observed that when the ship is in Case I using Case I 
optimum parameters it gave a better cost value than Case III 


parameters did, end of the 600 seconds simulation. BUG Vast 
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was also observed that if the final time of simulation is 
extended, the difference between the cost values decreases 
and if we continue to increase the final time, using Case 


III parameters in Case I gives a smaller cost value than 
Case I parameters do. 


These results show that the transient response part is 


important in finding the optimum control parameters. 
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VII. AN APPROACH TO AN ADAPTIVE AUTOPILOT 


p OE————————— (A ——- — — PH UO 


AS Seen from the previous chapter the optimal controller 


parameters change for changes in ship conditions such аз 


ship speed, loading etc. and also changes in environmental 
conditions such as sea state, encounter angle encounter 
frequency and depth of water. To maintain optimal steering 


performance automatically in the presence of changing condi- 
tions, it is necessary to design an adaptive system which is 
capable of self adjustment of the controller parameters to 


provide minimum added drag due to steering. 


The steering control system, if desired as an adaptive 
autopilot, would consist of four Subsystems as shown in 
Figures (gate) 

* Subsystem #1 would be a computer which will perform to 
find the optimal control parameters. "It ПОШ Кре ЫЕ 
information about ship steersngsecbaracterpg scr EN 
System state sensors such aS а gyrocompass, rudder 
angle potentiometer and speed log and it should feed 
the control parameter values to the controller. 

¢ Subsystem #2 would be the controller, which should be 
adjustable. TE gets the parameter values from 
subsystem #1 and sends the rudder command to subsystem 
#3 which steers the ship. 

* Subsystem #3 is the plant which includes the system 
State sensors and ship steering devices such as servo 
motors, hydraulic pumps, rudder, steering gear etc. 

* Subsystem #4 is a manual control option for safety 
rules. If manual steering is needed it cuts the 
connection between controller and steering devices and 


gives the control to the helmsman.In case of computer 
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failure it cuts the connection between computer апа 
controller and sends to the controller parameter values 
which are chosen by the watch officer. It contains 
look up tables which specify control parameters as 


functions of steering and environmental conditions. 


MANUAL 
CONTROL 
GROUP 


Comme TE Ik 


CONTROLLER 





Figure 7.1 Adaptive Control Scheme. 


\ 


Success of the this adaptive autopilot system depends on 
the computer program as well as the accuracy of the system 
sensors. The computer program which was used in this 
research may be used on board, but the present program mini- 
mization subroutine needs a lot of computation time and it 
also needs starting guesses for parameters which are 
different for every condition. If computation time is 


reduced to a reasonable time and starting values are made 


65 


available in proper limits, the modified function minimiza- 
tion program would be considered as an adequate program for 
on board purposes. The work on reducing the computation 


time is presented сл pet NM 


In the future, ships could have better measurement of 
navigation than can be provided by conventional equipment on 
board. For example, the U.S Navy is involved in a program 
to build a system which is called NAVSTAR/GLOBAL POSITION 
SYSTEM (GPS). The system will provide extremely accurate 
three-dimensional position and velocity information to users 
anywhere in the world. And also another system is called 
NAVY REMOTE OCEAN SENSING SYSTEM (NROSS) will be able to 
determine wind velocities over th world's oceans with an 
accuracy sufficient to determine ocean surface waves. Using 
such valid information the watch officer can use look up 
tables and insert them into the computer, so system opera- 
tion will be very close to the minimum cost value and the 
function minimization program will) .accomplish исе 
EuUDI tes Gap re ly Detailed information about GPS and NROSS 
can be found in [Ref. 10, 11, 12, 13]. 
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МЕРТ, CONCLUSTONS SANDS RECOMMENDABIONS FOR FUTURE STUDY 


A. 


CONCLUSIONS 


The conclusions resulting from this research of the 


Mariner class ship fuel consumption as it relate to steering 


might be listed as follows. 


A steering control system would minimize propulsion 
losses due to steering and maintain desired heading 
with reasonable heading error in every ship and enviro- 
mental condition. 

The study shows that the best model to describe the 
dynamics of the ship is the Taylor's series expansion, 
which allows both linear and nonlinear terms in the 
ship's equation of motions. Also the third order ship 
Nomoto model is reasonable to use instead of the ship's 
equation of motions. It involves both the sway and yaw 
equations. 

ТЕ is believed that the COS E ncCtUlomn which 5 
presented in Chapter 6 and is commonly used by many 
researchers is an adequate function for on board use. 
It has variables such as heading error and rudder angle 
which can be easily measured on board, and a weighting 
Factor (LAMBDA) which is also easy to calculate but 
depends on ship conditions such as ship speed. Results 
in chapter 7 show that accuracy of the LAMBDA value 
does not make significant changes in the controller. 
ies study two different types of controller were 
tried which have been called controller 'A' and 'B'. 
The structure of these controllers is shown in Chapter 
И Controller 'A' was determined to be a best struc- 
ture, because in some cases the adaptive calculations 


for controller 'B' did not converge. 
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• Ап adaptive controller which minimizes propulsion 
losses due to steering is needed when ship characteris- 
tics and environmental conditions change. 

* For performance in fuel saving, an adaptive controller 
may Бе better than the existing Universal  Gyropilot 
(UGP). 

e Since, equations of motion of Surface ships differ only 
in the numerical values of the coefficients, the simu- 
lation programs used for the Mariner class ship would 
be useable for other type of ships, knowing their hull 
characteristics.: The studies made tor Eshe SI m 


containership are examples.  [Ref. 6, 9, 7, 8]. 


B. RECOMMENDATIONS FOR FUTURE STUDY 


This research does not cover all ship and enviromentel 
conditions, so to get a better understanding about optimal 
controller parameters, it is recommended that the methods be 
applied to find the controller parameters for an expanded 


range of operating conditions. 


This thesis investigated only course keeping with 
emphasis on minimizing rudder and yawing activity to reduce 
fuel consumption. If a track following controls in autos 
matic control for replenishment at sea is desired, a 
different cost function might be needed. The nature of the 


cost function for such applications should be studied. 


Irregular seas case were not considered in this study. 
For future work it is necessary to study the effect of 
irregular seas on the controller parameters and on the cost 
function. ТЕ is believed that compering the regular sea 
results with irregular sea results would give better under- 


standing about ship's steering characteristics. 
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The transient response has a big affect on finding the 
optimum parameter values. For future studies it may be 
better to increase final time so that the effect of the 


transient response would not be significant. 


Recent studies оп roll stabilization shows that using 
ти ег stabilization is successful in reducing roll, 
[Ref. 14, 15]. Adding the roll equation into the ship model 
ama Getermining a proper cost function for minimum roll 
motion a controller could be designed fon roll 


свава Па 2 аетоти. 
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АРРЕМОТХ А 
THESIS FORTRAN 


$ JOB 


aa OQ 


REA 1, БО, О А о 
REAL*8 X,XDOT Y YDOT, Ü UDOT . Y , YDOT, YAV, R,RDOT 
REAL*8 TIME ERINE, X X2, X Х6,Х7,Х8 
REAL*8 YO,Yl,Y уа 74, 157 n^ 12 та“ 
ВЕА #8 МО, МТ, ND? МЗ, NB. NE NZ NS 
REAL*8 C1,C2,C3 Ch КЕ ЕЛЕ» 
REAL*8 RO,DELT,$ 45 Ul K П ГОРШ РО 
REAL*8 DYAÀWE,YÀWÉ ДИС? TSR ISE.: TDLFF , LAMDA 
REAL*8 51,5 А. 
REAL*8 MASS,IZ,.XG,YVDOT,NVDOT 
REAL*8 YR,YRDOT,NR,NRDOT,FX,FY,MZ 
REALX8 RX.RY КО ТХ, ТУ ТО, НА WE 
REAL*8 RYR,RXI,RYl,.MZR,MZI 
INITIAL CONDITIONS. FOR INTEGRATION 
SIMULATION END TIME IN SECONDS 
ETIME-210. 
COUN 
INITIALIZE THE COST FUNCTION 
SR=0.0 
TDTEF 020 
DA=2.91 
GAIN COEFFICIENTS TO BE OPTIMIZED 
217-450 
Pl-9.26 
X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 
ОКО 
DADOS 
YDOT-0.0 
U,UDOT,V.VDOT ARE FIX COORDINATES ON SHIP 
F2-0.0 
B3 OO 
ү-0.0 
UDOT-0.0 
VDOT-0.0 
YAW-0.0 
YAWE-0.0 
R-0.0 
RDOT-0.0 
ORDERED SPEED IN  FEET/SEC 
15.*1,689 FT/SEC=15 KNOTS 
01-15,%1,689 
AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
D = RUDDER ANGLE 
IUS 00572292 
Во. 
Ee pos 
L3j=L*L*L 
L4=L*L3 
L5=L*L4 
L6-L*L5 
FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
MOMENTS IN Z 
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м2-0 | 
RXR--.37126D5 
RXI-.68406D5 

ВУВ=-. 3998306 
RYI-.2447D6 

MZR-.296D8 
MZIz-.17546D8 i 
RX= (RXR**2+RXI**2 )* 
ВИ о рико 5 
RZ= (MZR**9+MZ1%*2 93 
TX=DATAN2 (RXI,RXR 
TY-DATAN2(RYI,RYI 
TZ-DATAN2(MZI,MZR 
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та 
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ПОЕНЕ 158+ 15К 
СО ТО 200 
400 CONTINUE 
STOP 
END 


ENTRY 
END 
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АРРЕМОТХ В 
DATA FOR SEA STATE PROGRAM 


The sea state program which is explained in [Ref. 4], 
needs information about the ship to calculate the ship's 
added mass, added inertia and seaway disturbance forces and 
moments. Information about the Mariner was gathered from 
[Ref. 27, 28] and presented here in the form which the sea 
State program needs. The current line is drawn on next page 


to show the Location of the values in cne о ај 
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APPENDIX C 
PROGRAM TO CALCULATE OPTIMAL GAINS 


Algorithm used here is showed in Figure ( C.1 ). 


UL ws D Е сша тате 
dH Controller of 


\ 
` 


+ motion 


Pune) lag 
Minimization 
SU pire Urine 


T DIFF 


Figure C.1 Algorithm to Calculate Optimal Gains. 


The disturbance forces (FX,FY) and moment (MZ) 


regular seas were aplied in program as a cosine wave. 


Precudurewtor des thice us: 
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for 
The 


* The real and imaginary value of force or moment is read 
from the sea state program output depending on ship 
Speed, encounter frequency and encounter angle. 

° These values are converted to magnitude and phase 


values. 


e Using the formula below, forces and moment are created 


Ma added into the ship s equations of motion. 
Force or Moment - WA*MAGNITUDE*cos(WE*TIME + PHASE) 


Where: 
WE = Encounter frequency (radian/second) 


WA = Significant wave height (ft) 


The correspondence between sea state and Significant wave 
Бп s = indicated 1п Table XXI [Ref. 29]. During this 
research the values that are presented in Table XXII were 


used as significant wave height (WA). 


TABLE XXI 
Sea State vs Range for Wave Height 


Sea State Range for wave height 

(Beaufort scale) (Feet ) 
0 0.0 
ili 0.32 
2 0.65-0.98 
3 1.96-3.28 
4 Dog Од 
5 6.56-8.20 
6 9.84-13.1 
7 1.18.2 
8 18.2-24.6 
9 25 0-32.9 
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TABLE XI 
Sea State vs Wave Height 





ea state Wave height (ft) 
5 5-30 
6 10.0 
/ 15.228 
8 2050 
9 2590 


The program is set up to calculate the optimal gains for 
controller A. It сап be modified to obtain optimal cains esac 
the rest of the controllers. After obtaining the optimal 
gains the program must be modified to do a simulation. The 


program has sufficient comments for appropriate changes. 


This program can be modified to obtain the Nomoto models. 


It is referenced in Chapter 4. The following need to be 
changed. 
C GAIN COEFFICIENTS TO BE OPTIMIZED 

K=XX(1) 

Zl-XX(2) 

Р1-ХХ(3) 

Р2-ХХ(4) 


C ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL - YAW COMMAND) 
c FOR EQUATIONS OF MOTION. 
D=YAW - YAWC 
ERROR SIGNAL TO DRIVE RUDDER (YAW COMMAND - YAW ACTUAL) 
C FOR NOMOTO 3RD ORDER MODEL. 
D2=YAWC- YAW2 
р01-(р2 - 01У)/РІ 
DQ2=((K* (2: DO) Om) ви 
C INTEGRATION 
Q1=Q1+DQ1*DELT 
Q2-Q2-*DQ2*DELT 
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YAW2-YAW2-Q2*DELT 
C COST FUNCTION 
TDIFF=TDIFF+ (YAW-YAW2)**2 


PROGRAM TO CALCULATE OPTIMAL GAINS FOR CONTROLLER 


//TANSAN JOB (1789,0356), RESEARCH' , CLASS-J 
/ | *MAIN ORG-NPGVMI.1789 
ЕЕЕ Ө ЕСЕ, ,IMSL- DP,REGION-1024K 
PEST SIN DD 
IN ORDER TO PERFORM SIMULATION,ONLY WHEN GAINS HAVE.BEEN 
" ,OBTAINED CHANGE XS(*) TORX (=) AND DELETE  XU(*),AND 


XL( 
DIMENSION XS(3),XU(3),XL(3) 
Х5(1)=0.6 
Х5(2)-41.58 
3)-11.33 
ST IS THE SSIaAPTING GUESS | 
C XL(I) IS THE LOWER LIMIT FOR THE I'TH VARIABLE 
ЕО „19 ТНЕ UPPER LIMIT FOR THE I'TH VARIABLE 
XU(1)=0.9 
RL (2 )=30. 
О) = 506 
XL(3)26. 
Х0(3)=15. 
С A DESCRIPTION OF THE FOLLOWING PARAMETERS 
C IS DISCUSSED IN BOXPLX 
6-2. (13. 
NTA= 1000 
NPR=0 
МАУ-0 
ЈЕ 
С THE FOLLOWING STATEMENT MUST BE CHANGED TO 
C CALL PLANT(XX 
C IF ONLY SIMULATION IS WANTED 
CALL BOXPLX(NV,NAV,NPR,NTA,R,XS,IP,XU,XL,YMN, IER) 
WRITE (6,25 
а FORMAT(1X,' OPTIMAL GAINS',/) 
РО 30 1=1.3 
30 WRITE(6 40)I SD | 
40 FORMA (J X 2. ПЕТА л) 


END 
SUBROUTINE TEN b 

C SUBROUTINE ктк SIMULATES ITHE SHIP 
COTON TDIFF 


REM ST L2.Lh3.L4,.L 


L6 
REAL XDOT,Y, TDT’ 5 UDOT i , VDOT YAW, В. ВОРОТ 
REAL TINE ETIME X1,X2,.X Х6,Х7,Х8 
REAL*8 YO,Yl,Y Кай 74, 98% 3% 12. 78 


МО, МТ, ИВ} МЗ, NA NS N6 N7 


DEDE 001.2263. C4À, Fl A ROS 
RO Сере? 5 ‚б Ul K F 
REAL*8 DYÀWE ,YÀWÈ VANC? 188 ISE? ЕВА LAMDA 
REAL: NDS p. 
REAL*8 MASS,.1Z,XG,YVDOT,NVDOT ,YR,YRDOT 


О 

NR, Naot, ХЕ» ,FY,MZ,RXR,RYR,RXI,RYI,MZR,MZI 
REAL" | I] TV T2 WA МЕ 

DIMENSION XX 


С INITIAL SION XS FOR INTEGRATION 
[E SINPEATTONSCEND TIME IN SECONDS 


J 
t] 
> 
# s , , , # 
Ыс оосо 


Sel 


С OO Qe 


ETIME=600. 


TIME=0.0 
ICOUNT= 1 
INITIALIZE THE COST FUNCTION 
повео О 
ISR-0.0 
TDIFF=0.0 
LAMDA=2.91 
GAIN COEFFICIENTS TO BE OPTIMIZED 
Ер t 
Pl=XX(3 
X, XDÖT, p OT ARE FIX COORDINATES ON EARTH 
EUN 
XDOT=9. 0 
ШИ а ARE FIX COORDINATES ON SHIP 
F2 QO 9 
F3=0.0 
WE о 
UDOT=0.0 
VDOT=0.0 
YAW=0.0 
YAWE=0.0 
EN 
RDOT=0.0 
ORDERED SPEED IN  FEET/SEC 
oil 689 ЕТ/5ЕС-15 KNOTS 
AT STEADY STATE ACTUAL SPEED (U) - COMMAND SPEED (UC) 
- RUDDER ANGLE 
D-0.0/57.296 
ео. 
ВЕР 
р 
ЕЕ. 
ВОВ 
L6=L*L5 
FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
MOMENTS IN Z 
EXT 
FY-0. 
М7=0. 
ВХВ=. 3337403 
ВХТ=. 1606403 
ВУВ=-. 1725415 
ВУТ=. 143504 
MZR-.30976D7 
Z1-.84672D6 " 
RX=(RXR**2+RXI**2)**,5 
RY-(RYR**2*RYI*2j*.5 
RZ-(MZR/*2*«MZI*2)7*.5 
TX-DATAN2(RXI,RXR 
TY-DATAN2(RYI,RYI 
TZ-DATAN2(MZI,MZR 
SIGNIFICA T WÀVE HEIGHT:(SEA STATE 2) 


W 

PNS EE FREQUENCY : (WHEN ENCOUNTER ANGLE IS 00) 
ADDED MASS AND ADDED INERTIA TERMS: 

MASS=. +07 


° ° — 
how Ovo nf 
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INTEGRATION 
=U+UDOT*DELT 


= TIME*DELT. 
Т-ІСООМТ»1 
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TIME,TDI 
500 FORMA .F9.2 


x. F20. 10) 


N 
DELETE ALL THE FOLLOWING SUBROUTINE SG: У ИО АНИНА 
AND NOT OPTIMIZATION IS WANTED 


SUBROUTINE BOXPLX (CATEGORY HO) 

PURPOSE 
BOXPLE IS A SUBROUTINE USED TO SOLVE THE PROBLEM OF 
LOCATING A MINIMUM (OR MAXIMUM) OF AN ARBITRARY OBJECT- 


IVE FUNCTION SUBJECT TO ARBITRARY EXPLICIT AND/OR 
IMPLICIT CONSTRAINTS BY THE СОМРШЕ Ж ОЕШ НИИ pro 
EXPLICIT CONSTRAINTS ARE DEFINED AS UPPER AND LOWER 
BOUNDS ON THE INDEPENDENT VARIABLES IMPLICIT CONSTRAINTS 
MAY BE ARBITRARY FUNCTION OF THE VARIE ESS ТВ 
CTION SUBPROGRAM TO EVALUATE THE OBJECTIVE FUNCTION AND 
IMPLICTIZGGNSIRAINIS., КЕЗБЕСТЕ Е, MUST PE SUPPLIED 

BY THE USER (SEE EXAMPLE BELOW BOXPLX ALSO HAS tHE 
OPTION TOSBEEBEORM INTECERSBISS IN WHERE THE VALUES 
OF THE INDEPENDENT VARIABLES 0 RESTRICTED ШОШ ТЕСЕК». 
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USAGE 
CALL BOXPLX (NV,NAV,NPR,NTA,R,XS,IP,XU,XL,YMN,IER) 


DESCRIET TON OF PARAMETERS 


NV AN_INTEGER INPUT DEFINING THE NUMBER OF INDEPENDENT 
КАКЕ РЕ ЖОЕ tne OBJECTIVE FUNCTION TO BE MINIMIZED. 
ЕЕ МУ + МАУ TS PRESENTLY 50. MAXIMIM NV IS 
ЕТ ВНЕ P TIMITS J ST ве EXCEEDED, PUNCH A SOURCE 
I S SH MES USUAL MANNER, AND CHANGE THE DIMENSION 


NAV AN INTEGER INPUT DEFINING THE NUMBER OF AUXILIARY var 
таре. BUIESUSERGWISHESCSTO DEFINE FOR HIS OWN CONVENIENCE. 
I O JP НЕ MAY WISH TO DEFINE THE VALUE OF EACH IMPLICI 
CONSTRAINT FUNCTION AS AN AUXILIARY VARIABLE. IF THIS 

IS DONE, THE OPTIONAL OUTPUT FEATURE OF BOXPLX CAN BE 
USED TO OBSERVE THE VALUES OF THOSE CONSTRAINTS AS THE 
SEBNEBISNEDROCGRESSES. S AUXILIARY VARIABLES, IF USED 

SI (POD ВЕ EVALUATED IN FUNCTION КЕ (DEFINED BELOW). 
М ВЕ ZERO. 


NPR INPUT INTEGER CONTROLLING THE FREQUENCY OF OUTPUT 
desired for ода nostic purposes. 

ТЕ ЕВ .LE. | s WILL BE 

PRODUCED BY ВОХР . OTHERWISE, THE CURRENT COMPLEX OF 
K= NV VERTICES AND THEIR CENTROID WILO BE OUTPUT AMEL 
EAC CHANPR PERMISSIBLE TRIALS. THE NUMBER OF TOTAL TRIAL 
NUMBER OF FEASIBLE TRIALS, NUMBER OF FUNCTION EVALUATIONS 
AND NUMBER OF IMPLICIT CONSTRAINT EVALUATIONS ARE IN- 
CLUDED IN THE OUTPUT. 

ADDITIONALLY (WHEN NPR „СТ. 0) THE SAME INFORMATION 
WILL BE OUTPUT: 


Weer THE INITIAL POINT IS NOT FEASIBL 

AE ШГЕК сети COMPLETE COMPLEX IS  CENERATED 

ШО ТЕ А OBEASBBLESVERTEX CANNOT BE FOUND AT SOME о ЕТАР, 

4) ТЕ ТНЕ ОВЈЕСТТУЕ VALUE OF A VERTEX CANNOT BE MADE 
NO-LONGER-WORST. 

3) IF THE LIMIIT ON TRIALS кт Т ne REACHED AND 

ODE AREN THE OBJECTIVE aas BEEN UNCHANGED FOR 

| ео INDICATING LOCAL MINIMUM HAS BEEN 


IF THE USER WISHES TO TRACE THE PROGRESS OF A лен, 
ДЕСИСТСЕ ОЕ МРК 25, 50 OR 100 IS RECOMMENDED 


NTA I INPUT OF LIMIT ON THE NUMBER OF TRIALS 
allowed ене calculation. 

ТЕ ТНЕ USER BNPEDSUNTA .EE.—0. A default 

VALUE OF 2000 IS USED. WHEN THIS LIMIT IS REACHED 
Copper Oe RETURNS TO THE CALLING PROGRAM WITH THE BEST 
ATTAINED OBJECTIVE FUNCTION VALUE IN YMN, AND THE BEST 
fT Ae РО SOLUTION POINT IN XS. 


R A REAL NUMBER INPUT TO DEFINE THE FIRST RANDOM NUMBER 
USED IN DEVELOPING THE INITIAL COMPLEX OF 2*NV VERTICIES. 
(0. fen. ee ГТ. 1.2 IF R TŠ NOT WITHIN THESE BOUNDS, 

T WILL BE REPLACED BY 1./3. 


Х5 INPUT REAL ARRAY "¿usa AT LEAST МУТМАМ. 

them irsi ny must cont 

FEASIBLE ORIGIN FOR STARTING THE CAL- 

GULATION.. THE LAST NAV NEED NOT BE INITIALIZED. UPON 
m TROM BOXPIX, THE FIRST NV ELEMENTS OF THE ARRAY 
CONTAIN THE COORDINATES OF THE MINIMUM OBJECTIVE 
function, AND THE REMAINING NAV (NAV .GE. 63 СОМТАТМ ТНе 
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оао ооо оосо ОООО ОООО ОООО ооа 
C) 


values of THE CORRESPONDING AUXILIARY VARIABLES. 


ШЕ INTEGER INPUT FOR OPTIONAL INTEGER PROGRAMMING. 
if ip-1, THE VALUES OF THE INDEPENDENT We PEP Vee 
be replaced WITH INTEGER VALUES (STILL STORED AS REAL*4). 


XU A REAL ARRAY DIMENSIONED AT LEAST NV INPUTTING THE 
upper BOUND ON EACH INDEPENDENT VARIABLE CEU EXPLIGCII 
conSTRAINT). INPUT VALUES ARE SLIGHT ALTERED BY BOXPLX. 


XL PREAL ARRAY DIMENSIONED AT LEAST NV INPUTTING THE 
lower bound each а A 

VARIABLE (ЕАСН EXPLICIT CONS tz Е) 

NOTE: FOR BOTH XU AND XL NOOSE. REASONABLE 

VALUES IF NONE ARE GIVE NOT VALUES WHICH AR 

magnitudes ABOVE OR BELOW THE EXPECTED SOLUTION. 

input values are SLIGHTLY ALTERED BY BOXPLX. 


YMN THIS OUTPUT IS THE VALUE (REAL*4) OF THE OBJECTIVE 
funcTION,CORRESPONDING TO THE SOLUTION POTNRR I TrOT ies 


ТЕК INTEGER ERROR RETURN. ТӨ БЕ INTEREOGA Pap eer on 
return FROM BOXPLX. IER WILL BE ONE OF ІНЕ FOREOWING: 


=-] CANNOT FIND FEASIBLE VERTEX OR FEASIBLE CENTROID 
AT THE START OR A RESTART (SEE 'METHOD' BELOW). 

=0 FUNCTION VALUE UNCHANGED FOR 'N' TRIALS. (WHERE 
N=6*NV+10 THIS IS ТНЕ МОЮ, КЕТШЕМ РАКАМЕТЕЕБЕ 

l CANNOT DEVELOP FEASIBLE VERTEX, 
2 CANNOT DEVELOP A NO-LONGER-WORST VERTEX, 

3 LIMIT ON TRIALS REACHED. NTA EXCEEDED) 
OTE: VALID RESULTS MAY BE RETURNED IN ANY ÓF THE 

ABOVE CASES. 


EXAMPLE OF USAGE 
THIS EXAMPLE MINIMIZES IHE.OBJECPRIVESRUNGITONSSHODLEIISR 
the T a RD) FE(X). THERE ARE TWO INDEPENDENT 
VarIABLE S TU. M s ще TWO IMPLICIT CONSTRAINT 
function ү! al QUE EVALUATED AS AUXILIARY 
variables (see S. NAL FUNCTION КЕ ХО, 
DIMENSION XS(4),XU(2),XL(2) 


SIS NDS EL 


ди пиг 


= 
Ex) 
КУ а га 
са 
H 


Ба 

= 
РА РА 
а ССП 

и 
N 
Wu WMO S—— c... 
H 
H 


l H XH H I I 

H 
OW OCOrHnocHOrm 
е. Ге е Ге е 
OO OO {л 


оон 
СО. 


Х (REUS NAV ;NPR, NTA, К, XS „ТР, XU,XL,YMN,IER) 
ИИ) RE PÓI e TS ТЕ Ана = је 
x 

E #6 CTION VALUE 18 T! FEP то, 
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e 


ОООО ООООООО 020000 Q0O000000000 00000 000010000 O0OOcaoodO3909003000000000000000 ОООО 


со 


FUNCTION ЕНЕ 
EVALUATE CONSTRA S. S 
1 ОТОК ЗЕТ KE-1 
constraint is violated. 
em) X(4) 


| 


| 

2 

ет. 

=з oo 
И Т 
4) „СЕ. 


Е ВТР МО ШИРЕТЕТ CONSTRAINT 
TEENY IMPLICIT 


ра 
Ni 
ии 


71320515X 
D TOURS 5% 3) ИМ О СОТО 1 
з Ел 


RETURN 


FUNCTION ЕЕ/Х) 
DIMENSION X( 


THIS IS THE OBJECTIVE FUNCTION 
TETUER 
END 


METHOD 


BHESDSMPLEX MEPFHOD IS e EXTENSION AND ADAPTION OF 

the simple method of linear рвота: 

STARTING WITH ANY ONE feasible point in n-dimension 
COMBEEX OF 2*N vertices is constmucted by 

SELECTING RANDOM POINTS WITHIN THE feasible 

КЕШЕ в 15 PURPOSE N COORDINATES ARE FIRST 

RANDOMLY CHOSEN WITHIN THE SPACE BOUNDED BY EXPLICIT CON- 

STRAINTS. THIS DEFINES A TRIAL INITIAL VERTEX. 

it is then checked for poss А s ер 

GE ЛЫР ТСТ CONSTRATNTS. ТЕ ore are violated, 

THE TRIAL INITIAL VERTEX IS DISPLACED half of its 

DISTANCE EBOBMOTHECCENTROID OF PREVIOUSLY SELECTED initial 

VERTICES: ВВЕЛА 115 DISPLACEMENT. PROCESS 15 

КЕНЕ БЕР ЧЛ НЕ VERTEX HAS BECOME FEASIBLE. IF THIS 

fails to happen after 5*n+10 P EAD. 

THE SOLUTION IS ABANDoned. after each vertéx is added 

TO SGHE. COMEEE X THE CURRENT centroid is checked for 


02) /(46.76538)) 


FEASIBILITY. İF IT IS INFEASIBLE, the last trail 
VERTEX IS ABANDONED AND AN EFFORT 10 GENERATE an alter- 
ATIVE TRIAL VERTEX IS MADE. 5*N-«10 VERTICES ARE 


ABANDONED CONSECUTIVELY , THE SOLUTION ШОТ ЕЕЙТМАТЕР. 


ТЕ АМ с e up Да ТОНЕВ, fake SAC 

computa 

TAE E INSTRUCTI NS. FIND. THE CURRENT WORST vertex that 
THE VERTEX МЕН QE LARGEST CORRESPONDING value КОТ 


ЕДЕ OBJECTIVE FUNCTIO AND REPLACE THAT VERTEX BY 
ITS OVERSE REFLECTION THROUGH THE CENTROID OF ALL OTHER 
vertices. if the vertex to be 


REPLACED I CONSIDERED AS A VECTOR IN n-space, 

ITST TOVE TST TEFLECTION IS OPPOSITE IN DIRECTION, IN 
CREASBEBNTNOLENGTH BY THE FACTOR 1.3, AND COLLINEAR WITH 
THE REPLACED VERTEX AND CENTROID OF ALL OTHER VERTICES. 


WHEN AN OVER- REFLECTION IS NOT FEASIBLE OR REMAINS 
worst, it is Considered not permissible 

AND IS DISPLACED HALFWAY TOW RD the centroid. 

AFTER FOUR SUCH ATTEMPTS ARE MADE UNSUCCESSFULLY 

EVERY FIFTH ATTEMPT IS MADE BY О THE OFFENDING 


87 


ОООООООООО ОООООООООДООООО ОООООООО ООООО ОООООСКЛО ООООООООООООО ООДОДООО СООООО 


vertex through the present best 
VERTEX. INSTEAD OF THROUGH THE CENtroid. if 5*n+10 
DISPLACEMENTS AND OVER-REFLECTIONS OCCUR without a 
SUCCESSFUL M RESULT, THE CURRENT TEESI 
pd: IS TAKEN AS AN INITIAL FEASIBLE POINT FOR A 
start run of the complete process. 
RESTARTING 15 ALSO UNDERTAKEN when 6*nv-*10 consecutive 
TRIALS HAVE BEEN MADE WITH NO SIGNIFicant 5 їп the 
VALUE OF THE OBJECTIVE FUNCTION. IN ALL case 
RESTARTING TS INHIBITED LF THE LAST RESTART DID NOT 
Е SIGNIFICANT IMPROVEMENT IN THE MINIMUM 
attaine 


IT IS RECOMMENDED THAT THE USER БЕТ ЕВЕ НОВ РОВ 
FURTHER USEFUL INFORMATION. 11 SHOUED ВЕЕТ ЕНЕ 
ALGORITHM DEFINED THERE НА5 ВЕЕМ АБТЕКЕО ОИЕ 
CONSTRAINED MINIMUM, RATHER THAN THE MAXIMUM. 


REMARKS 


THE INTEGER PROGRAMMING OPTION WAS ADDED TO THIS PROGRAM 

AS SUE RID IN REFERENCE. (2). А MIXED 

inte continuous variable version of а 

WOUL "B EASY TO CREATE BY DEclarin o be an array 

OF NV CONTROL VARIABLES WHERE IP ur POT: indicate 

e THE I-TH VARIABLE IS TO РБЕ Г М D EO Eo D 
ALUE > EACH STATEMENT OF THE мерс 

WOULD THEN NEED TO BE ALTERED TH Е {ТР Е etc. 
WHERE THE SUBSCRIPT IS ERN LLY, 

XU AND E ARE ALTERED TOTBETAN EESTPON ` WITHIN 

actua Da 

DECLARED BY THE USER. THIS ADJUSTMENT IS NOT MADE 

WHEN IP-1. 


NOTE: МО NON-LINEAR PROGRAMMING ALGORITHM CAN GUARANTEE 
that the answer found is the globa 

MINIMUM, RATHER 1 JUST À local minimum. however, 
ACCORDING TO REF.2, THE COMPLEX method has an advantage 
IN THAT IT TENDS тд FIND THE GLOBAL minimum more 
FREQUENTLY THAN MANY OTHER NON-LINEAR PROGRAM- 
MING ALGORITHMS. 


Т ДРВНА EN p пох VARIABLE FEATURE 
can a 

PROBLEMS CONTAINING EQUALITY CONstraints. any equality 
CONSTRAINT TIMPLIES Tar A GIVEN VARiable is not truly 
INDEPENDENT. THERE IN GENERAL, ONE variable 

INVOLVED IN AN EQUALITY CONSTRAINT CAN BE RENUMBERED from 
THE SET OF NV INDEPENDENT МАКТАВГЕРФИЛРРЕГИТТИШПЕ СО 
OF NAV AUXILIARY VARIABLES THIS USUALLY INVOLVES 
renumbering THE INDEPENDENT VARLABLES OF THE GIVEN 


le 
SUBROUTINES AND FUNCTIONS REQUIRED 


SUBROUTINE 'BOUT' AND FUNCTION 'FBV' ARE INTEGRAL 
parts of THE BOXPLX PACKAGE. 


Duo AU ие MUST BE SUPPLIED BY ТИЕ ШС ПНВ. 
is used to evaluate the implicit 
STRAINTS. SET KE=0 AT THE oR traint Of hunetd o 

ATHEN EVALUATE THE IMPLICIT s s s S. in the example 
м THE FIRST CONSTRAINT X(3) must be within the 
RANGE (0. ВЕ X(3 NEN É SECOND constraint x(4) 

MUST DO. IF eft ER CONSTRAINT IS pot үт лат 
THESE BOUNDS CONTR OL IS TRANSFERRED TO S ERENT 
AND KE IS SET TO "1" AND CONTROL IS RETURNED TO BOXBLX. 
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AOC C) COMO TO 


THE SECOND QEUNCTION THE ,USER MUST PROVIDE EVALUATES THE 
ща ective n 
LLED E AS SHOWN. IN THE PUE above. and fe 
TO THE VALUE OF THE OBJECTIVE function 
ОРС TO CURRENT VALUES OF THE NV INDEPENDENT 
VARIABLES IN ARRAY 'X' 
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ЕЕ ВОХРЫХ (NV. NAV, NPR,NIZ,RZ,XS,1P,BU,BL,YMN,IER) 


DIMENSION V(50,50), FUN(50), SUM(25), CEN(25), XS(NV) 
l,bu(nv),bl(NV) 


КУ = 5 

EP = 1.E-6 

NTA = 2000 

IF (NI2. О) МГА = №2 

IF (R.LE.0..OR.R.GE.1.) R-1./3. 


NVT 2 NV*NÀV 
TOTAL VARS, EXPLICIT PLUS IMPLICIT 
CURRENT TRIAL NO. 
CURRENT NO. OF PERMISSIBLE TRIALS 
TFCURRENT NO. OF TIMES F HAS BEEN ALMOST UNCHANGED 
CHECK FEASIBILITY OF START POINT 
-] ХУ 
E LE.VT) GO TO 1 
Г 


МРТ = 


r 


Cohn» 


I)-AMAXI(EP,EP^ABS 


C103 O3 H4 C) «Á HA SÁ HA HALO EA HL SUC 
СО "dpi nHmHiBH"odmm/''biHO 


Aa EP: FABS (BL(T))} 


I 
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GOO AA NAA 


aa 


Фе С) Cima 


NCE = l 
NUMBER OF CONSTRAINT EVALUATIONS 


KE(V( r 1 ВО Со TOSS 
ТЕ ( E11 3 95949718 
WRITE S 
GOTO 12 
5 NFE = 
NUMBER ОЕ VERTIGES (K) = 2 TIMES NO. OF VARIABLES. 
K = 2*NV 
NUMBER OF DISPLACEMENTS ALLOWED. 
NLIM = 5*NV+10 
NUMBER OF CONSECUTIVE TRIALS WITH UNCHANGED FE TO 
terminate. 
NCT = NLIM+NV 
ALPHA = 1.3 
FK = K 
FKM = FK-1. 
BETA = ALPHA+1. 
INSURE . ae RANDOM NUMBER GENERATOR IS ODD. 
(MOD(IOR, 2). ЕО. 0) IQR=IQR+101 
SET UP INITIAL VERTICES 
POOD = РУС, ша) 
ТД Е ща 
FUNOLD = FUN(1) 
DO 15 1° 2s 
FI - 


T 
TEN ји 
LIMT - 0 
7 ЕЕ ВТ Е 


END CALCULATION IF FEASIBLE CENTROID CANNOT BE FOUND. 
IF (LIMT.GE.NLIM) GO TO 11 


DO 8 J=1,NV 
“TOR „198 GENERATOR (RANDU) 
(IQR E IQR = IQR+2147483647+1 


Ë X 
R X = 656613E-9 
у s IE * ROX* (BU 3 ER 
I .EQ.1 (J.L )=AINE 6)? 
8 CON INUE 
DO_10 L= 1,NLIM 
NCE СЕ+1 
ТЕ E I))2EQ-.0) еони ние 
vr 2 Е v CEN 
+ 
ТЕ (IP. НЕЙ У Т+.5) 
и aly 
9 C NTT UE 
10 CONTINUE 
Ii ШЕ {ат LE-0) COTO 


NEUE iut NPT, NFE, NCE, NVS NVT V TEONE CENTI) 
GO RTO 48 


12 
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с) OQ 


ССС С 


С) 


OOA ССС) 


Ce) ©. OO 


OOO) 


13 DO 14 J-l,N 
и ИЗ) ие. 


TRY TO ASSURE FEASIBLE CENTROID FOR STARTING. 
IF (KE(CE 0) GO TO 60 
ВЫ E Sones j 


0) 
60 МЕЕ = 
15 FUN(T) = иа, 1)) 
END OF LOOP SETTING oF INITIAL COMPLEX. 
IF (NPR.L 


CALE BOUT Е 03 69 МРТ, ТИРЕ. МСЕ , МУ, МУТ ,У,К,ЕОМ , СЕМ, О) 
FIND THE WORST VERTEX, THE 'J'TH. 
ШО 16 Т2 К 
IF (FUN(J).GE.FUN(I)) GO TO 16 
16 CONTINUE 
BASIC LOOP. ELIMINATE EACH WORST VERTEX IN TURN. 
it must become NO LONGER WORST NOT MERELY IMPROVED. 
find иска to vertex, THE 'JN' TH. O 
ЕТ 
Е ВО. Г) JN = 2 
DO 18 i= 1, 
IF и EQ: Це соо 
IF ( FUN(JN).GE.FUN(L)) GO TO 18 
18 CONTINUE 
LIMT = NUMBER OF MOVES DURING THIS TRIAL TOWARD THE 
centroid DUE TO FUNCTION VALUE. 
DET = I 
COMPUTE CENTROID AND OVER REFLECT WORST VERTEX. 


Но 13 Ут ; J), 


SUM(I) a “ү 

СЕМ(Т) = k [FKM | 
ЕТА - ALPHA*VT 

F (IP.EQ. CEN р AINT(VT+.5) 


INSURE THE EXPLICIT CONSTRAINTS ARE OBSERVED. 
19 V(I,J) s» AMAX1(AMINI(VT,BU(I)),BL(I)) 


Ne = NT+1 
CHECK FOR IMPLICIT CONSTRAINT VIOLATION. 
20 DO 25 N=1,NLIM 
ТЕ“ EEC). ЕО 0) 5001026 
EVERY 'KV'TH TIME .QVER-REFLECT THE OFFENDING VERTEX 


through the BEST EX. 
IF [MOD QD(N EV] NE. 40) GÓ TO 22 


CALL 


р0 21 151, 
УТ = ВЕТА: AI, M)-ALPHA*V (I,J) 


2] 


OOQ a 


OOO A 


бс 


OO 


~“ OO OO 


СІСІСІСІ OQA 


AFO , AMX , FUNTRY ае мою, NTFS,N 


21 С. БО МАХ] ТАМ AMIN DEPO), BL DD 


СОТО 
CONSTRAINT VIOLATION: MOVE NEW POINT TOWARD CENTROID. 
22 DO 23 I: 
а M 
ТЕ (тр. Е V AINT(VT+.5) 
23 өзіні ЧЕ 


24 NT = NT+l 
25 CONTINUE 


IER = 1 
CANNOT GET FEASIBLE VERTEX BY MOVING ДД CENTROID, 
ORTEY OVER- Е THRU THE EEST S vEP [E 

TF (NPP .LE 0) GO TO 2 


52 
ATE BU ue NPT, NFE,NCE,NV,NVT,V,K, FUN, CEN, J) 


JEER EEE VERTEX FOUND, EVALUATE THE OBJECTIVESFUOUNCTTONS 
FUNTRY = FE(V(1, Jy) 
TEST IO SEE IE FUNGTTON р HAS NOT CHANGED. 
АЕО = A FUNTRY-FUN o 
Х = АМАХ1 (АВ$ (ЕР* PUNO D), EP) 
ACTIVATE THE FOLLOWING TWO STATEMENTS FOR DIAGNOSTIC 
purposes only. 
(67220 
99 FORMAT (X I El 
ТЕ (АЕО: С AMKS ес eee 27 
NTF NEES E 
1E Mu ЕТ. МСТ СО ТО 5 
F (NPR, ТЕ: 006 GO TO 42 
WRI Me ates 
CALL B Ur (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN, CEN, O) 
27 МТЕЗ = 0 


IS THE NEW VERTEX NO LONGER WORST? 
28 IF (FUNTRY.LT.FUN(JN)) GO TO 34 


TRIAL VERTEX IS STILL WORST: АВТ ТАРО СРО. 
EVERY TH TIME OVER- REFLECT THE OFFENDING VERTEX 
through th the BEST VER 


IF (MODI LIMT D. NE. 0) CORTO 
CALL FBV (K, ТУ” М) 


рО I=1 
УТ = ВЕТА* М) - -ALPHA? NC, J 
Ws p. EQ. 1 (19D; : ТУ 

29 XL (AMIN CV BUT JBE EI) 
GO TOMON 


30 DO 31 IiljNV 
Т T1 NI а ат, 5) 


22 


оо 00 о 


OO 


OQO OOQ C C) ОЮ ОСС 


ОО ССС 


ІРІСІ 


31 CONTINUE 
Вт. Г. МТМ) СО 1033 


САММОТ МАКЕ ӨТ = "J'TH VERTEX NO LONGER WORST BY 
displacin 


Trin SG TER = 2 OR" BY OVERSREFEECTING THRU THE BEST VERTEX. 


(N PR LE. OD GOTTO 42 
WRI E oor" 2) TAN 
CALL BO BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN, J) 


33 NT = кті 
GO TO 20 


PUGCESS: WE ME 7 LEP EACEMENT ЕРЕ ЕКЕ НЕ 
34 CS N OL) = PONT 
= FUNTRY 
NPI ESSET Sp 


EVERY 100'TH PERMISSIBLE TRIAL, RECOMPUTE CENTROID 
summation to AVOID CREEPING ERROR. 
IF (MOD (NPT. 100). ЫШ ОСО ТОУ 


DO 36 I=1,NV 


SUM(I) = 0. 
ШО Ет 
35 т. = SUM(I)*V(I,N) 
SUM(I)/FK 
36 CON ТИЕ (т) 
MENO 
GO TO 39 
37 DO 38 I-1,NV 
38 SUM(I) = SUM(I)*V(I,J) 
В 


39 IF DT LE. 0) 50 
IF (MOD(NPT,NPR). NE. 09 СОЛТО 0 


ООР КО Т ЫРТ,МЕЕ,МСЕ,МУ,МУТ,У,К,ЕОМ,СЕМ EC) 
HAS THE MAX. NUMBER OF TRIALS BEEN REACHED WITHOUT 
convergence? iF NOT, GO TO NEW TRIAL. 

40 IF (NT.GE.NTA) GO po j 

= көт VERTEX NOW BECOMES WORST. 
GO TO 1 

ТЕК 3 

ТЕ (МРВ. GT.0) WRITE (6,54) 


COLLECTOR POINT FOR ALL ENDINGS. 


41 


1 CANNOT DEVELOP FEASIBLE VERTEX. IER l 
Д CANNOT DEVELOP A NO-LONGER-WORST VERTEX. IER = 2 
3 FUNCTION VALUE UNCHANGED FOR K TRIALS. IER = Q 
4 LIMIT ON TRIALS REACHED. ТЕК = 3 
5 CANNOT FIND FEASIBLE VERTEX AT START. TER -I 
42 CONTINUE 


FIND BEST VERTEX. 
LL FBV (К.Р FUN M) 
I? ius GELS) GO TO 44 
RESTART IF THIS SOLUTION IS SIGNIFICANTLY BETTER THAN 


99 


ОО 


s IP (NPR LE ОБЕ d E. IS THE PIR iia 
(е QU (M, YMN FUN (M 
43 IE (Fu н) Е. MN) GO T | 
АВ5 N(M)- TM ).ЪЕ.АМАХ1(ЕР,ЕР УР БЕРСЕ? 
CINE TIT UD TRY UNLESS LIMIT ON ТЕТЕ. 


44 a = FUN(M 
FUN(1) = UN (M) 


l H H H 
< 


&5 V(I 
46 


GO TO 48 


ШЕ? LED) G GO TO 6 
UN CM ше NCE, NV, NVT. V, K, FUN V CIE TDR) 


47 


wri 
ae 
p 
Uri 
„си 


m E ка 56) ЕЏ 


49 FORMAT UM AND DIRECTION OF OUTLYING 
lvariable at star 
50 FORMAT SOHOIMPLICI CONSTRAINT VIOLATED Sar 


51 FORMAT ('0C OCANNOT. IND FEASIBLE',I4,'TH VERTEX OR 
Conrro-c St 

52 FORMAT (LOHOAT | TRIA I4,54H CANNOT FIND FEASIBLE 

lvertex which no 
SLONGER WORST. 14. 15X,'RESTART FROM BEST VERTEX. ") 
ара HAS BEEN ALMOST UNCHANGED 
га115 

54 PORMAT’ (27HOLIMIT ON TRIALS EXCEEDED. 

55 FORMAT OBEST VERTE? 1 IS NO.',I3,' OLD MÍN WAS',El15.7, 
56 FORMAT ("OMIN EA FUNCTION IS ',E15.7) 


SUBROUTINE FBV LIN 
DIMENSION FUN(50) 
М = 1 
DO 1 I= 
та (FUNC) Е PUN GI) C Os R. 
1 CONTINUE 
RETURN 
D 


МУ ТЕ) 


4 
É) F IX GUN 1, МУ 
disi N( І), (Y (3 (СРИИ р = ) 
| WRITE 6) ud 1) ,J=NVP,NVT) 

IF (IK.NE.0) GO TO 2 
WRITE кс D (COLI PN) 
0) с 


3 
пи (6, Ра) ет) I=1;NV) 


94 


3 WRITE (6,9) IK,(C(I),I-1,NV) 
RETURN 


4 FORMAT ('ONO. TOTAL и = T5. 4X. 
т. feasible trails = 15, 
25 NO: FUNCTION EVALUATIONS - 115 E 
3 по con traint шн 
A FUNCTION VALUE' ЕВЕ ENT VARIABLES/ 
5dependENT OR IMPLICIT До нае NTS' 
5 depend ТН ‚Е18, 7, 2X ,7E14. 2x, VEU NEN 
6 FORMAT 21X,.7E 2) 
7 FORMAT TOHOCENTRO ШЕН ОЛЕ Л. 7 (215 7Е14. 2), 
8 FORMAT ('0 BEST VERTEX ,7Х,7Е14 71215, ТЕЛЕ ИІ) 
9 ОСЕШ ОТО LESS QX то ок, / (21X, 


DIMENSION X(3 

COMMON TDIF 

CALL PLANT(X) 

FE-TDIFF 

RETURN 

END 

FUNCTION KE n) 

DIMENSION X(3 

KE=0 

RETURN 

END | 
/|[GO.SYSIN DD * 


THIS CARD SHOULD BE USED ONLY REGULAR SEA CASE 
//GO.FT12F001 DD DISP-SHR,DSN-MSS.S2160.A341 


FUNCTION ЊЕ) 
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